function [ L ] = Constructmatrix( Nx, Ny, Ux, Uy, dx, dy ,h ,g )
%Construction of the Spatial Matrix

dx1=1/dx;
dx2 = 1/(dx*dx);

dy1= 1/dy;
dy2 = 1/(dy*dy);


% ---- Construct the spatial Discretization Matrix------%

%------------- Construct the Hyperbolic Term--------------%

A = generate_hyperbolic( Nx, Ny, Ux*0.5*dx1, Uy*0.5*dy1); 
% A = zeros(Nx*Ny,Nx*Ny);
% % Inner Boundary points
% for j=2:Nx-1
%     for i=2:Ny-1
%         index = (j-1)*Ny+i;
%         A(index, index-1)= -Uy*0.5*dy1;
%         A(index,index+1) = Uy*0.5*dy1;
%         A(index, index+Ny) = Ux*0.5*dx1;
%         A(index,index-Ny) = -Ux*0.5*dx1 ;
%     end
% end
% 
% % East and West Boundary
% for i=2:Ny-1
%     j=1;
%     index = (j-1)*Ny+i;
%     A(index, index-1)= -Uy*0.5*dy1;
%     A(index,index+1) = Uy*0.5*dy1;
%     A(index, index+Ny) = Ux*0.5*dx1;
%     A(index,(Nx-1)*Ny+i) = -Ux*0.5*dx1 ;
%     
%     j=Nx;
%     index = (j-1)*Ny+i;
%     A(index, index-1)= -Uy*0.5*dy1;
%     A(index,index+1) = Uy*0.5*dy1;
%     A(index, i) = Ux*0.5*dx1;
%     A(index,index-Ny) = -Ux*0.5*dx1 ;
% end
% 
% % SW Corner point
% i=1;
% j=1;
% index = (j-1)*Ny+i;
% A(index, index+Ny-1)= -Uy*0.5*dy1;
% A(index,index+1) = Uy*0.5*dy1;
% A(index, index+Ny) = Ux*0.5*dx1;
% A(index,(Nx-1)*Ny+i) = -Ux*0.5*dx1 ;
% 
% % SE Corner Point
% j=Nx;
% index = (j-1)*Ny+i;
% A(index, index+Ny-1)= -Uy*0.5*dy1;
% A(index,index+1) = Uy*0.5*dy1;
% A(index, i) = Ux*0.5*dx1;
% A(index,index-Ny) = -Ux*0.5*dx1;
% 
% 
% 
% % North and South Boundary
% for j=2:Nx-1
%     i=1;
%     index = (j-1)*Ny+i;
%     A(index, index+Ny-1)= -Uy*0.5*dy1;
%     A(index,index+1) = Uy*0.5*dy1;
%     A(index, index+Ny) = Ux*0.5*dx1;
%     A(index,index-Ny) = -Ux*0.5*dx1 ;
%     
%     i=Ny;
%     index = (j-1)*Ny+i;
%     A(index, index-1)= -Uy*0.5*dy1;
%     A(index,index-Ny+1) = Uy*0.5*dy1;
%     A(index, index+Ny) = Ux*0.5*dx1;
%     A(index,index-Ny) = -Ux*0.5*dx1 ;
% end
% 
% % NW Corner Point
% i=Ny;
% j=1;
% index = (j-1)*Ny+i;
% A(index, index-1)= -Uy*0.5*dy1;
% A(index,index-Ny+1) = Uy*0.5*dy1;
% A(index, index+Ny) = Ux*0.5*dx1;
% A(index,(Nx-1)*Ny+i) = -Ux*0.5*dx1 ;
% 
% % NE Corner Point
% 
% i=Ny;
% j=Nx;
% index = (j-1)*Ny+i;
% A(index, index-1)= -Uy*0.5*dy1;
% A(index,index-Ny+1) = Uy*0.5*dy1;
% A(index, i) = Ux*0.5*dx1;
% A(index,index-Ny) = -Ux*0.5*dx1 ;
%------------- Construction of the Hyperbolic Matrix Finished--------------%

%------------- Compile the Matrices into Global Matrix L--------------%

L=zeros(2*Nx*Ny,2*Nx*Ny);
L(1:Nx*Ny,1:Nx*Ny)=A;
L(Nx*Ny+1:2*Nx*Ny,Nx*Ny+1:2*Nx*Ny)=A;

%------------- Construction of the Elliptic Term --------------%
% A = zeros(Nx*Ny,Nx*Ny);
% 
% % Populate the diagonal term
% for j=1:Nx;
%     for i=1:Ny
%         index = (j-1)*Ny+i;
%         A(index,index) = -2*h*(dx2+dy2);
%     end
% end
% 
% % Inner Boundary points
% for j=2:Nx-1
%     for i=2:Ny-1
%         index = (j-1)*Ny+i;
%         A(index, index-1)= h*dy2;
%         A(index,index+1) = h*dy2;
%         A(index, index+Ny) = h*dx2;
%         A(index,index-Ny) = h*dx2 ;
%     end
% end
% 
% 
% 
% % East and West Boundary
% for i=2:Ny-1
%     j=1;
%     index = (j-1)*Ny+i;
%     A(index, index-1)= h*dy2;
%     A(index,index+1) = h*dy2;
%     A(index, index+Ny) = h*dx2;
%     A(index,(Nx-1)*Ny+i) =h*dx2 ;
%     
%     j=Nx;
%     index = (j-1)*Ny+i;
%     A(index, index-1)= h*dy2;
%     A(index,index+1) = h*dy2;
%     A(index, i) = h*dx2;
%     A(index,index-Ny) = h*dx2 ;
% end
% 
% % SW Corner point
% i=1;
% j=1;
% index = (j-1)*Ny+i;
% A(index, index+Ny-1)= h*dy2;
% A(index,index+1) = h*dy2;
% A(index, index+Ny) = h*dx2;
% A(index,(Nx-1)*Ny+i) =h*dx2 ;
% 
% % SE Corner Point
% j=Nx;
% index = (j-1)*Ny+i;
% A(index, index+Ny-1)= h*dy2;
% A(index,index+1) = h*dy2;
% A(index, i) = h*dx2;
% A(index,index-Ny) = h*dx2;
% 
% 
% 
% % North and South Boundary
% for j=2:Nx-1
%     i=1;
%     index = (j-1)*Ny+i;
%     A(index, index+Ny-1)= h*dy2;
%     A(index,index+1) = h*dy2;
%     A(index, index+Ny) = h*dx2;
%     A(index,index-Ny) = h*dx2 ;
%     
%     i=Ny;
%     index = (j-1)*Ny+i;
%     A(index, index-1)= h*dy2;
%     A(index,index-Ny+1) = h*dy2;
%     A(index, index+Ny) = h*dx2;
%     A(index,index-Ny) = h*dx2;
% end
% 
% % NW Corner Point
% i=Ny;
% j=1;
% index = (j-1)*Ny+i;
% A(index, index-1)= h*dy2;
% A(index,index-Ny+1) = h*dy2;
% A(index, index+Ny) = h*dx2;
% A(index,(Nx-1)*Ny+i) = h*dx2;
% 
% % NE Corner Point
% 
% i=Ny;
% j=Nx;
% index = (j-1)*Ny+i;
% A(index, index-1)= h*dy2;
% A(index,index-Ny+1) = h*dy2;
% A(index, i) = h*dx2;
% A(index,index-Ny) = h*dx2 ;

A = generate_elliptic(Nx,Ny,h*dx2,h*dy2);
%------------- Compile the Matrices into Global Matrix L--------------%
L(1:Nx*Ny,Nx*Ny+1:2*Nx*Ny)=A;

for k=1:Nx*Ny
    L(Nx*Ny+k,k)=g;
end
end

